Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.
In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv
Model code is in ./models/production/modelProduction_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/modelProduction_OB_makeFile.R
Functions for this qmd file in ./models/qmdProduction_OB_functionsToSource.R
There are 3 models here: All cohorts:
model 1, with quadratic weight effect on phi By cohort:
model 2, with linear weight effect on phi model 3, with quadratic weight effect on phi
19.0.1 How many ageInSamples to include?
Code
all <-tar_read(cdWB_electro_target)#str(all)table(all %>%filter(river =="wb obear")|> dplyr::select(ageInSamples))
19.1 All Cohorts, modelNum 1: phi(i,t,c) * g(i,t,c), p(i,t) model
19.1.1 Model code
Code
# all cohorts have the same code for a given model. show model code for one of them hereload('./models/production/runsOut/production_OB_model_1_mostRecent.RData') # WILL NEED TO CHANGE TO MOST RECENT#load('./models/production/runsOut/production_OB_model_2_2025-02-28 19.RData')d$modelCode
{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]],
sd = 2)
weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = weightSD)
for (t in (first[i]):(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]/90
weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1],
sd = weightSD)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] <- grIntT[t] + grBeta[1, cohort[i]]
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ dnorm(0, sd = 0.5)
}
for (c in 1:nCohorts) {
grBeta[1, c] ~ dnorm(0, sd = 1/0.667)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i,
t] * weight[i, t] + phiBeta[2, i, t] * weight[i,
t] * weight[i, t] + phiBetaCohort[cohort[i]]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
for (b in 1:2) {
phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
}
}
}
for (c in 1:nCohorts) {
phiBetaCohort[c] ~ dnorm(0, sd = 1/0.667)
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
phiIntTOut[t] <- ilogit(phiIntT[t])
phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
}
for (i in 1:N) {
for (t in 1:(T - 1)) {
weightZ01[i, t] <- weightIfAvailable(weight[i, t],
sdWeight, meanWeight, z[i, t], available[i, t])
z01[i, t] <- aliveIfAvailable(z[i, t], available[i,
t])
}
}
for (t in 1:(T - 1)) {
sumWeightZ01[t] <- sum(weightZ01[1:N, t])
sumZ01[t] <- sum(z01[1:N, t])
}
for (cohortLoop in 1:nCohorts) {
for (t in 1:(T - 1)) {
sumWeightZ01ByCohort[cohortLoop, t] <- sumVecForCohort(weightZ01[1:N,
t], cohortLoop, cohort[1:N])
sumZ01ByCohort[cohortLoop, t] <- sumVecForCohort(z01[1:N,
t], cohortLoop, cohort[1:N])
}
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}
sumWeightZ01ByCohortS <-MCMCsummary(object = d$mcmc, params =c("sumWeightZ01ByCohort" ), round =3) %>%rownames_to_column(var ="var")|>mutate(cohort0 = var |>str_extract("\\[(\\d+),") %>%# Changed to match number before commastr_remove_all("\\[|,") %>%# Remove [ and commaas.numeric(),ais = var |>str_extract(",\\s*(\\d+)\\]") %>%str_remove_all(",|\\s|\\]") %>%as.numeric(),cohort = cohort0 + cohorts[1] )sumZ01ByCohortS <-MCMCsummary(object = d$mcmc, params =c("sumZ01ByCohort" ), round =3) %>%rownames_to_column(var ="var") |>mutate(cohort0 = var |>str_extract("\\[(\\d+),") %>%# Changed to match number before commastr_remove_all("\\[|,") %>%# Remove [ and commaas.numeric(),ais = var |>str_extract(",\\s*(\\d+)\\]") %>%str_remove_all(",|\\s|\\]") %>%as.numeric(),cohort = cohort0 + cohorts[1] )
Code
ggplot(sumZ01ByCohortS, aes(ais, mean, group = cohort, color = cohort)) +geom_line(aes(linetype =factor(cohort))) +labs(#title = "sumZ01 by cohort and age in samples",x ="Age in samples",y ="Estimated number of fish alive")
Code
ggplot(sumWeightZ01ByCohortS, aes(ais, mean, group = cohort, color = cohort)) +geom_line(aes(linetype =factor(cohort))) +labs(#title = "sumZ01 by cohort and age in samples",x ="Age in samples",y ="Estimated sum of fish mass")
19.2 By Cohort, modelNum 2: phi(i,t) * g(i,t), p(i,t) model
Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.
Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)
library(targets)library(nimble)# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html# # # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAssource('./models/production/modelProduction_OB_functionsToSource.R')source('./models/production/qmdProduction_OB_functionsToSource.R')modelNum <-2# phi * growth
19.2.2 Model code
Code
# all cohorts have the same code for a given model. show model code for one of them hereload('./models/production/runsOut/production_OB_model_2_2002_mostRecent.RData')d$modelCode
{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]],
sd = 2)
weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = weightSD)
for (t in (first[i]):(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]/90
weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1],
sd = weightSD)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ dnorm(0, sd = 1000)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i,
t] * weight[i, t]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
for (b in 1:1) {
phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
}
}
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
phiIntTOut[t] <- ilogit(phiIntT[t])
phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
}
for (i in 1:N) {
for (t in 1:T) {
weightZ01[i, t] <- weightIfAvailable(weight[i, t],
sdWeight, meanWeight, z[i, t], available[i, t])
z01[i, t] <- aliveIfAvailable(z[i, t], available[i,
t])
}
}
for (t in 1:T) {
sumWeightZ01[t] <- sum(weightZ01[1:N, t])
sumZ01[t] <- sum(z01[1:N, t])
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}
#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same ehojs_define(summary_OB_OJS_all_2 = sAll_2)#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))
Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
19.2.4.3 Plot survival
Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Combine scohort$summaryMod2_tt_growth data into one data frame
Code
sParamsAll0_2 <-tibble(.rows =0) |>add_column(!!!s2002_2$summaryMod3_tt_growth[0,]) # need to change to mod2# Fill tibblefor (cohort in cohorts) { summary_obj_2 <-get(paste0("s", cohort, "_2"))$summaryMod3_tt_growth # need to change to mod2 summary_obj_2$cohort <- cohort sParamsAll0_2 <-bind_rows(sParamsAll0_2, summary_obj_2)}# remove brackets for filteringsParamsAll0_2$varNoIndex <- sParamsAll0_2$var |>str_remove("\\[.*\\]")# for numeric orderingsParamsAll1_2 <- sParamsAll0_2 |>filter(varNoIndex !="phiBetaT") |>mutate(varNumeric1 = var |>str_extract("\\[(\\d+)\\]") %>%str_remove_all("\\[|\\]") %>%as.numeric() )sParamsAll2_2 <- sParamsAll0_2 |>filter(varNoIndex =="phiBetaT") |>mutate(varNumeric1 = var |>str_extract("\\[(\\d+),") %>%# Changed to match number before commastr_remove_all("\\[|,") %>%# Remove [ and commaas.numeric(),varNumeric2 = var |>str_extract(",\\s*(\\d+)\\]") %>%str_remove_all(",|\\s|\\]") %>%as.numeric() )sParamsAll_2 =bind_rows(sParamsAll1_2, sParamsAll2_2)
19.2.5.1 Combined cohorts
19.2.5.2 p
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.2.5.3 phi
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="phiIntTOut"), aes(varNumeric1, mean, group = cohort, color =factor(cohort))) +geom_point() +geom_line() +xlab("AgeInSamples")
19.2.5.4 grIntT
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="grIntT"), aes(varNumeric1, mean, group = cohort, color =factor(cohort))) +geom_point() +geom_line() +xlab("AgeInSamples")
19.2.5.5 grIntT, filtered
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="grIntT", mean >-5, mean <5), aes(varNumeric1, mean, group = cohort, color =factor(cohort))) +geom_point() +geom_line() +xlab("AgeInSamples")
19.2.5.6 phiBetaT
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="phiBetaT", mean >-5, mean <5), aes(varNumeric2, mean, group = cohort, color =factor(cohort))) +geom_point() +geom_line() +xlab("AgeInSamples") +facet_wrap((~varNumeric1))
19.2.5.7 sumWeightZ01
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color =factor(cohort))) +geom_point() +geom_line() +xlab("AgeInSamples")
19.2.5.8 sumZ01
Code
ggplot(sParamsAll_2 |>filter(varNoIndex =="sumZ01"), aes(varNumeric1, mean, group = cohort, color =factor(cohort))) +geom_point() +geom_line() +xlab("AgeInSamples")
19.2.6 Size-dependent survival
Need to get rid of squared term for cohort-specific models
19.3 By Cohort, modelNum 3: phi(i,t) * g(i,t)^2, p(i,t) model
Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.
Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)
library(targets)library(nimble)# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html# # # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs#source('./models/production/modelProduction_OB_functionsToSource.R')#source('./models/production/qmdProduction_OB_functionsToSource.R')modelNum <-3# phi * growth
19.3.2 Model code
Code
# all cohorts have the same code for a given model. show model code for one of them hereload('./models/production/runsOut/production_OB_model_3_2002_mostRecent.RData')d$modelCode
{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]],
sd = 2)
weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = weightSD)
for (t in (first[i]):(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]/90
weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1],
sd = weightSD)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ dnorm(0, sd = 1000)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i,
t] * weight[i, t] + phiBeta[2, i, t] * weight[i,
t] * weight[i, t]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
for (b in 1:2) {
phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
}
}
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
phiIntTOut[t] <- ilogit(phiIntT[t])
phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
}
for (i in 1:N) {
for (t in 1:T) {
weightZ01[i, t] <- weightIfAvailable(weight[i, t],
sdWeight, meanWeight, z[i, t], available[i, t])
z01[i, t] <- aliveIfAvailable(z[i, t], available[i,
t])
}
}
for (t in 1:T) {
sumWeightZ01[t] <- sum(weightZ01[1:N, t])
sumZ01[t] <- sum(z01[1:N, t])
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}
#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same ehojs_define(summary_OB_OJS_all_3 = sAll_3)#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))
Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
19.3.4.3 Plot survival
Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Combine scohort$summaryMod2_tt_growth data into one data frame
Code
sParamsAll0_3 <-tibble(.rows =0) |>add_column(!!!s2002_3$summaryMod3_tt_growth[0,]) # Fill tibblefor (cohort in cohorts) { summary_obj_3 <-get(paste0("s", cohort, "_3"))$summaryMod3_tt_growth summary_obj_3$cohort <- cohort sParamsAll0_3 <-bind_rows(sParamsAll0_3, summary_obj_3)}# remove brackets for filteringsParamsAll0_3$varNoIndex <- sParamsAll0_3$var |>str_remove("\\[.*\\]")# for numeric orderingsParamsAll1_3 <- sParamsAll0_3 |>filter(varNoIndex !="phiBetaT") |>mutate(varNumeric1 = var |>str_extract("\\[(\\d+)\\]") %>%str_remove_all("\\[|\\]") %>%as.numeric() )sParamsAll2_3 <- sParamsAll0_3 |>filter(varNoIndex =="phiBetaT") |>mutate(varNumeric1 = var |>str_extract("\\[(\\d+),") %>%# Changed to match number before commastr_remove_all("\\[|,") %>%# Remove [ and commaas.numeric(),varNumeric2 = var |>str_extract(",\\s*(\\d+)\\]") %>%str_remove_all(",|\\s|\\]") %>%as.numeric() )sParamsAll_3 =bind_rows(sParamsAll1_3, sParamsAll2_3)
19.3.5.1 Combined cohorts
19.3.5.2 p
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.3.5.3 phi
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.3.5.4 grIntT
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="grIntT"), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.3.5.5 grIntT, filtered
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="grIntT", mean >-5, mean <5), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.3.5.6 phiBetaT
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="phiBetaT", mean >-5, mean <5), aes(varNumeric2, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples") +facet_wrap((~varNumeric1))
19.3.5.7 sumWeightZ01
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.3.5.8 sumZ01
Code
ggplot(sParamsAll_3 |>filter(varNoIndex =="sumZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +geom_point() +geom_line() +xlab("AgeInSamples")
19.3.6 Size-dependent survival
Need to get rid of squared term for cohort-specific models